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SUMMARY 

A  two-step  approach  to  finite  element  ordering  is  introduced.  The  scheme  involves  ordering  of 
the  finite  elements  first,  based  on  their  adjacency,  followed  by  a  local  numbering  of  the  nodal 
variables.  The  ordering  of  the  elements  is  performed  by  the  Cuthill-McKee  algorithm.  This 
approach  takes  into  consideration  the  underlying  structure  of  the  finite  element  mesh,  and  may  be 
regarded  as  a  "natural"  finite  element  ordering  scheme.  The  experimental  results  show  that  this 
two-step  scheme  is  more  efficient  than  the  reverse  Cuthill-McKee  algorithm  applied  directly  to  the 
nodes,  in  terms  of  both  execution  time  and  the  number  of  fill-in  entries,  particularly  when  higher 
order  finite  elements  are  used.  In  addition  to  its  efficiency,  the  two-step  approach  increases 
modularity  and  flexibility  in  finite  element  programs,  and  possesses  potential  application  to  a 
number  of  finite  element  solution  methods. 

INTRODUCTION 

The  use  of  finite  element  methods  typically  involves  solving  a  system  of  equilibrium  equations  : 
K  u  *  R  (1) 

where  u  and  R  are,  respectively,  the  displacement  and  the  loading  vectors.  K  is  the  global 

stiffness  matrix  which  is  often  symmetric,  positive  defimte,  and  populated  with  many  zeros.  In 

the  solution  process,  it  is  important  to  take  into  account  the  structure  of  the  matrix  K.  i.e.  the 

pattern  of  the  zero  and  nonzero  entries  so  that  the  computational  effort  is  minimized. 

The  solution  process  for  solving  Equation  1  can  be  divided  into  four  separable  tasks1*  : 

1.  Ordering  -  to  find  a  proper  permutation  matrix  P,  such  that  the  symmetric  permuted 
matrix  PKP1  has  a  desirable  structure. 

2.  Storage  allocation  -  to  determine  the  necessary  information  about  the  structure  of  the 
matrix  factor  of  PKP1  and  to  set  up  the  storage  scheme. 

3.  Numeric  factorization  -  to  decompose  the  permuted  matrix  PKP'  into  a  triple  matrix 
product  LDL1. 

4.  Solution  -  to  compute  the  solution  vector  u  successively  by  a  forward  substitution, 
z  *  L*'(PlR),  and  a  backward  substitution,  u  *  PtL^D^z). 

The  identification  of  these  four  separable  tasks  not  only  encourages  software  modularity,  but  also 
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facilitates  the  theoretical  study  of  sparse  matrices.  This  paper  focuses  on  the  first  task,  that  of 
ordering  of  a  system  of  equilibrium  equations  resulting  from  finite  element  discretization. 

When  the  matrix  K  is  decomposed  into  its  matrix  factors  L  and  D.  it  typically  suffers  some 
fill-in:  that  is.  the  filled  matrix.  F  *  L  +  L1,  has  nonzeros  in  positions  which  are  zero  in  K.  One 
objective  in  ordering  a  matrix  is  to  minimize  the  number  of  these  fill-in  entries.  By  doing  so. 
the  computing  cost  may  also  be  reduced  since  the  amount  of  computation  depends  on  the  number 
of  nonzeros  m  the  matrix  factor  L. 

A  finite  element  mesh  consists  of  a  collection  of  finite  elements  which  are  connected  at  their 
common  boundaries.  The  discretization  of  the  mesh  and  the  finite  elements  are  selected  so  that 
the  behavior  of  the  physical  structure  under  examination  can  be  properly  modelled.  The  nodal 

variables  are  then  defined  on  the  finite  element  mesh;  their  number  depends  on  the  type  of  finite 
elements  used.  Traditionally,  most  ordering  schemes  have  been  developed  for  ordering  directly  the 
nodal  variables  and.  hence,  the  equilibrium  equations  in  1.  Very  little  effort  has  been  reported  m 
developing  ordering  schemes  which  include  consideration  of  the  underlying  topological  structure  of 
the  finite  element  mesh. 

One  solution  scheme  specifically  designed  for  solving  a  system  of  equations  resulting  from  finite 
element  problems  is  the  frontal  solution  method14.  This  method  is  often  considered  as  the  most 
"natural”  solution  scheme  since  it  operates  directly  on  the  underlying  structure  of  the  finite 

element  mesh.  In  this  solution  scheme,  the  finite  elements  are  assembled  and  entered  into  the 
solution  one  at  a  time.  A  nodal  variable  is  eliminated  as  soon  as  all  the  neighboring  fmite 

elements  incident  to  it  are  available.  The  remaining  assembled  nodal  variables,  which  are  not  yet 
ready  for  elimination,  are  retained  in  the  "front",  or  "active  front". 

Another  solution  scheme  which  is  similar  to  the  frontal  solution  method  is  the  so-called 

generalized  element  method20.  In  this  method,  the  nodal  variables  of  an  active  front  are  treated 

as  a  superelement.  Multiple  fronts  are  allowed,  and  each  front  is  treated  as  a  superelement. 

This  method  has  recently  attracted  attention  from  numerical  analysts6.  The  merit  of  the  frontal 

solution  and  the  generalized  element  methods  is  that  both  can  be  extended  to  include 
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considerations  of  auxiliary  storage  in  a  simple  manner. 

Obviously,  the  efficiency  of  the  frontal  solution  and  the  generalized  element  methods  depends  on 
the  order  in  which  the  finite  elements  are  numbered.  The  ordermg  of  the  nodal  variables  is 
implicitly  imposed  by  the  order  in  which  the  finite  elements  are  processed.  To  the  authors’ 
knowledge,  there  have  not  been  any  studies  reported  for  ordering  the  finite  elements  directly.  In 
this  paper,  we  attempt  to  examine  the  feasibility  of  ordering  the  finite  elements. 

It  will  be  shown  that  ordering  the  finite  elements  has  considerable  computational  advantage  over 
ordering  the  nodal  variables,  but  that  ordering  the  finite  elements  alone  does  not  completely  define 
the  permutation  matrix  P.  It  is  therefore  necessary  to  include  into  the  ordermg  scheme  a  local 
ordering  strategy  for  numbering  the  nodal  variables  of  each  finite  element  so  that  the  number  of 
fill-in  entries  can  be  minimized.  This  "two-step"  approach,  which  includes  the  ordermg  of  the 
finite  elements  as  well  as  of  the  nodal  variables,  is  the  subject  of  this  paper. 

GRAPH  THEORY  NOTATION  AND  SPARSE  MATRIX  STUDY 
The  use  of  graph-theoretic  approaches  has  found  many  applications  m  sparse  matrix  study. 
Although  few  results  from  graph  theory  are  directly  applicable  to  the  study  of  sparse  matrices, 
the  graph  representation  is,  nonetheless,  a  powerful  tool  to  characterize  the  structure  of  a  sparse 
matrix. 

Symmetric  Matrices  and  Finite  Graphs 

A  finite  undirected  graph  G=(V,E)  consists  of  a  finite  set  V  of  n  elements,  called  nodes ,  and 

a  set  E  of  unordered  pairs  of  distinct  nodes  (v.v),  called  edges.  For  any  node  v  in  G,  the  set 

■  j 

of  nodes  adjacent  to  v,  adj(v),  is  defined  as 

adj(v)  =  {  v  (  V  |  (v.v  )  (  E  }. 

i  i 

The  number  of  nodes  in  adj(v),  denoted  by  ladj(v)!,  is  called  the  degree  of  v.  The  deficiency 
of  v,  def(v),  is  the  set  of  distinct  pairs  of  nodes  in  adj<v)  which  are  not  themselves  adjacent.  A 
graph  is  complete  if  every  pair  of  nodes  is  adjacent.  A  subgraph  G'-( V'.E’)  of  G  is  a  graph 
for  which  V’  c  V  and  E’  c  E.  A  complete  subgraph  is  called  a  clique.  A  section  graph 
G(V’)  is  a  -subgraph  G*(V',E(V’))  induced  by  a  node  set  V’,  where 
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E(V’)  =  {  (v,v)  (  E  !  v  ,v  <  V’} 

>  j  '  j 

If  a  set  of  nodes  V”  is  deleted  from  the  graph  G*<  V,E),  the  section  graph  G(V\V")  is  obtained 

from  G  by  removing  the  nodes  V"  together  with  their  incident  edges. 

A  (simple)  path  {v.,  v }  is  a  sequence  of  distinct  nodes  and  continuous  edges  leading  from 
v  to  v  such  that  there  are  no  repeating  edges.  A  graph  G  is  said  to  be  connected  if  there  is 
at  least  one  path  between  every  pair  of  distinct  nodes  in  G;  otherwise,  G  is  disconnected.  A 

connected  subgraph  is  called  a  component.  The  distance  d(v,v)  between  two  nodes  v  and  v 

'  j  1  j 

in  a  connected  graph  is  the  number  of  edges  in  the  shortest  path  joining  nodes  v  and  v.  The 

eccentricity  e(v)  of  a  node  v  is  then  given  by 
e(v)  *  mas  {  d(v,v)  !  v  «  V  }. 

i  i 

The  diameter  5(C)  of  the  graph  G=( V.E)  is  defined  as 
5(0  *  max  f  e(v),  v  <  V  }. 

A  node  v  is  said  to  be  a  peripheral  node  if  e(v)  =  3(G). 

For  a  graph  G»(V,E)  with  n  nodes,  an  ordering  (or  labelling)  of  V  is  a  bijection  (a  one-to- 

one  onto  mapping)  : 

*  :  {  1,  2 .  n  }  <->  V. 

The  ordered  graph  of  G  is  denoted  by  Gfl  *  <V,E.«).  The  integer,  ranging  from  1  to  n, 
assigned  to  a  node  by  the  ordering  is  called  the  number  or  label  of  that  node. 

The  notation  of  a  level  structure  is  commonly  used  to  describe  the  properties  of  many 

ordering  schemes.  A  level  structure  LS  of  a  graph  G  is  an  arrangement  of  the  nodes  of  G  into 

m  levels,  ordered  LS .  LS  ,  such  that  nodes  at  a  given  level  LS  are  connected  to  no  nodes 

1  m  e  i 

other  than  LS  ,  LS  and  LS  .  For  a  node  v  <  V,  there  corresponds  a  level  structure  rooted 

k-l  k  k-1  r 

at  v.  The  levels  of  the  rooted  level  structure  are  determined  by 

1.  assigning  LS^  =  {v},  and 

2.  for  k  >  1,  assigning  to  LS^  the  set  of  nodes  adjacent  to  nodes  in  LS^  (  which  have 
not  yet  been  included  in  any  previous  levels. 

Given  an  n  by  n  symmetric  matrix  K,  there  corresponds  to  it  a  finite  graph  G*(  V.E),  where  a 
node  v  «  V  denotes  the  i<h  row  (or  column)  of  the  matrix  K.  and  an  edge  (v,v>  t  E 

1  >  i 
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symbolizes  an  offdiagonal  nonzero  entry  K.  .  t*  K  ).  An  ordering  of  the  graph  G  of  matrix  K.  is 
equivalent  to  a  symmetric  permutation  PKP1,  where  P  is  an  n  by  n  permutation  matrix.  The 
unordered  graphs  (in  which  nodes  are  not  labelled)  of  K  and  PKP'  are  the  same,  but  the  node 
numbers  of  the  associated  ordered  graphs  are  different. 

Finite  Element  Mesh  and  Its  Solution  Graph 

A  finite  element  mesh  is  a  collection  of  finite  elements  m  which  the  adjacent  finite  elements 
are  joined  at  their  common  boundaries  or  vertices.  There  is  a  node  at  each  vertex  of  the  finite 
element  mesh.  For  finite  elements  of  higher  order,  where  higher  order  polynomial  interpolation 
functions  are  used,  there  may  also  be  nodes  lying  along  the  sides  or  on  the  faces  of  the  finite 
elements,  and/or  located  internally  within  the  finite  element  itself.  The  nodes  situated  at  the 
mterior  of  the  finite  element  are  referred  as  the  interna/  variables .  and  the  nodes  located  along 
the  boundaries  are  called  the  external  variables. 

Associated  with  each  node  is  a  set  of  variables  corresponding  to  the  degrees  of  freedom  defined 
at  that  node.  This  set  of  variables,  in  turn,  corresponds  to  a  submatrix  m  the  global  stiffness 
matrix,  which  is  typically  full.  During  assembly  and  solution,  this  submatrix  can  conveniently  be 
treated  as  a  single  entity.  For  simplicity,  the  set  of  variables  at  a  node  is  referred  as  the  nodal 
variable. 

The  finite  element  mesh  can  be  transformed  directly  into  the  finite  graph  G  representing  the 
structure  of  the  global  stiffness  matrix.  We  call  this  graph  a  finite  element  solution  graph. 
or,  simply,  a  solution  graph  (to  distinguish  it  from  the  finite  element  connectivity  graph  to  be 
introduced  below).  The  nodes  of  the  solution  graph  are  the  nodal  variables  defined  on  the  finite 
element  mesh.  The  edges  of  the  solution  graph  are  constructed  by  making  the  nodes  of  each 
finite  element  pairwise  adjacent,  since  a  nonzero  entry  K  in  the  global  stiffness  matrix  K.  implies 
that  the  i1*1  and  the  jlb  nodal  variables  share  at  least  one  incident  finite  element.  In  other  words. 


if  V’  denotes  the  set  of  nodal  variables  belonging  to  one  finite  element,  the  section  graph  G(V') 
of  the  solution  graph  G  is  a  clique:  the  union  of  the  cliques  over  all  elements  defines  the  edges 
in  the  solution  graph.  Examples  of  transforming  a  finite  element  mesh  into  its  associated  solution 
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graph  are  shown  in  Figure  1. 

Matrix  Factorization  and  Its  Graph  Theoretic  Model 

The  factorization  of  K  can  be  symbolically  modelled  using  a  graph-theoretic  approach.  This 

approach  is  particularly  helpful  in  understanding  how  the  fill-m  entries  are  created  during 

factorization. 

The  most  fundamental  scheme  of  matrix  factorization  is  one  analogous  to  Gaussian  elimination, 
summarized  in  Figure  2.  At  the  ilh  step  of  the  factorization,  the  ith  column  of  L  and  the  ilh 
diagonal  entry  of  D  are  computed,  and  the  matrix  K*1’  is  condensed  and  modified  to  K<l*n.  It  is 
in  the  condensation  that  the  fill-in  entries,  if  any,  are  created. 

Using  the  finite  graph  representation  of  a  sparse  matrix,  the  factorization  scheme  can  be 

modelled  graph-theoretically  as  a  node-elimination  process1*.  Upon  elimination  of  node  v,  the 

elimination  graph  G  is  obtained  from  C-(V.E)  by  : 

1.  deleting  node  v  and  its  incident  edges: 

2.  adding  auxiliary  (fill-m)  edges  such  that  all  adjacent  nodes  of  v  form  a  clique. 

That  is, 

G  *  ( V\v,  E(V\v)  U  def(v)). 

V 

For  an  ordered  graph  Ga,  the  elimination  is  a  recursive  process  defined  by 
P(G  >  *  {  G  1  G  G .  G  } 

a  o  “  1  n-l 

where  G  is  called  the  i,h  elimination  graph  obtained  bv  eliminating  v  from  G  The  fill-m 

1  I  l  —  l 

edees  created  dunne  the  elimination  of  v.  denoted  bv  f,  are  the  deftv  )  in  G  These  fill-m 

-  -  i  i  ii-l 

edges  correspond  to  the  new  nonzero  entries  introduced  into  the  condensed  matrix  K"*'1  at  the  i‘b 
step  of  the  factorization. 

In  a  f;nite  element  solution  graph.  the  new  clique  formed  by  the  elimination  of  a  node  can  be 
treated  as  a  new  finite  element,  called  a  generalized  element  or  a  superelement20.  The  node- 
elimination  process  on  a  finite  element  graph  can  thus  be  interpreted  as  a  series  of  transformation 
of  the  finite  elements  into  the  superelements.  This  elimination  model  also  suggests  an  interesting 
characteristic  related  to  higher  order  finite  elements.  Internal  nodes  of  a  higher  order  element  are 


i 
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connected  only  to  the  external  nodes  of  that  element.  Furthermore,  all  nodes  in  a  finite  element 
form  a  clique.  Hence,  if  an  internal  node  is  eliminated  before  any  of  the  external  nodes  of  the 
same  element,  there  is  no  fill-in  since  the  deficiency  of  the  internal  node  is  null. 

After  elimination,  the  set  of  fill-in  edges  is  given  by 
n-1 

4><G  )  =  U  f. 

*  »=1  ' 

The  filled  graph  G^,  which  represents  the  structure  of  the  filled  matrix.  F  *  L  +  L1.  is  defined 
by 

G  «  (V,  E  U  <t><G  )). 

F  a 

The  node-elimination  model  was  formally  introduced  by  Rose18.  Since  then,  efforts  have  been 
made  to  develop  efficient  algorithm*;  implementing  the  model19,  These  algorithms  have  been 
encoded  into  sparse  matrix  programs  such  as  YMSL'  and  Sparspak9.  The  major  application  of 
these  algorithms  is  in  the  second  task  presented  in  the  Introduction,  namely  to  set  up  the  data 
structure  for  stormg  the  numeric  entries  of  the  matrix  factors. 

A  node-addition  model  which  simulates  the  Cholesky  factorization  has  recently  been  developed 
by  the  authors15.  For  this  graph-theoretic  model,  nodes  are  added  onto  the  filled  giaph  G^,  rather 
than  being  eliminated  from  the  original  graph  Gfl.  The  structure  of  the  matrix  factor  L  is 
constructed  one  row  at  a  time.  When  computing  the  entries  of  a  particular  row  in  L.  the  model 
does  not  require  any  a  priori  information  for  the  rows  beyond  the  current  row.  Therefore,  this 
model  has  the  flexibility  that  the  tasks  of  labelling  a  nodal  variable,  determining  its  row  structure 
m  L,  and  computing  the  numeric  entries  of  the  row  can  all  be  performed  simultaneously. 

The  execution  times  for  the  numeric  factorization  as  well  as  the  symbolic  factorization  of  a 
given  matrix  K  depend  on  the  number  of  nonzeros  m  the  matrix  factor  L.  Therefore,  it  is 
worthwhile  to  develop  efficient  algorithms  to  order  the  matrix  K.  so  that  the  number  of  fill-in 
entries,  or  equivalently  the  number  of  nonzeros  m  L,  is  minimized.  In  the  next  section,  two  well 
known  ordering  algorithms,  namely  the  Cuthill-McKee  algorithm  and  the  reverse  Cuthill-McKee 
algorithm,  are  discussed. 
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The  Cuthill-McKee  and  the  Reverse  Cuthill- McKee  Ordering  Algorithms 

Let  K  be  an  n  by  n  symmetric  matrix.  For  the  ith  row  of  K,  define 

a  (K)  3  min  {  j  I  K  ^  0  } 

1  u 

and 

B  (K.)  3  l  *  a  <K). 

1  I 

The  number  a  (K)  denotes  the  column  subscript  of  the  first  nonzero  entry  in  the  1th  row  of 

matrix  K.  The  number  ^(K>  is  usually  referred  as  the  local  bandwidth  of  the  i*  row  of 

matrix  K.  and  is  equal  to  the  number  of  off-diagonal  entries  between  the  first  nonzero  of  the 

row  and  the  mam  diagonal.  The  bandwidth  and  the  proiile  of  matrix  K  are  defined  as 
b(K)  =  max  {  fi  (K),  i  3  1,  ...n  } 

and 

n 

/>< K>  3  I  fi  (K) 

i=l  ’ 

respectively.  The  nonzero  entries  of  the  filled  matrix  F  of  K  are  confined  within  the  local  band 
of  each  row.  Many  ordering  schemes  have  been  developed  to  minimize  the  bandwidth  or  the 
profile  of  K. 

One  ordering  scheme  commonly  used  with  finite  element  programs  to  minimize  the  bandwidth  is 
the  Cuthill-McKee  algorithm4.  Basically,  the  algorithm  is  a  breadth-first  technique  in  labelling 
the  nodes  of  a  graph.  Once  a  finite  element  mesh  is  tranformed  into  its  finite  element  solution 
graph,  the  algorithm  numbers  the  nodes  in  the  following  way. 

1.  Determine  a  starting  node  and  number  it  v^. 

2.  For  nodes  v,  i=l„..,a,  find  all  unlabelled  neiehborine  nodes  of  the  node  v  and 

i  w  i 

number  them  sequentially  in  increasing  order  of  degree. 

This  ordering  algorithm  is  essentially  the  same  as  that  for  generating  a  level  structure  rooted  at 
the  starting  node  v  . 

In  selecting  a  starting  node,  a  peripheral  node  is  preferred13.  The  idea  of  choosing  a  peripheral 
node  as  a  starting  node  is  to  generate  as  many  levels  as  possible  m  the  level  structure,  since  an 
increase  in  the  number  of  levels  tends  to  decrease  the  profile  of  the  corresponding  matrix. 
Unfortunately,  existing  algorithms  for  finding  a  peripheral  node  are  not  computationally  feasible. 
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A  compromise  that  has  been  suggested  is  to  choose  a  starting  node  with  high  eccentricity.  This 
node  is  generally  referred  as  a  pseudo- peripheral  node.  An  efficient  algorithm  to  find  a 
pseudo-peripheral  node  has  been  suggested  and  implemented  by  George  and  Liu10.  For  most 
practical  cases,  the  execution  time  for  finding  a  pseudo-peripheral  node  is  no  greater  than 
0(!E1),  i.e.,  the  order  of  the  number  of  edges  in  the  graph  G=(V,E). 

George  eLal.12  have  shown  that  if  linear  insertion  is  used  for  sorting  the  time  complexity  for 

the  second  step  of  the  Cuthill-McKee  algorithm  requires  ai  most 
(  4 1 E  |  ♦  2cm  1 E I  )  operations. 

where  !E!  is  the  number  of  edges  in  the  graph  G,  m  is  the  maximum  degree  of  any  node,  and  c 
is  some  constant. 

It  has  been  discovered  that  significant  improvements  in  minimizing  the  profile  can  be  achieved 
by  simply  reversing  the  node  ordering  obtained  from  the  Cuthill-McKee  algorithm*.  The  resulting 
algorithm  is  the  well  known  reverse  Cuthill-McKee  (RCM)  algorithm,  which  can  be  summarized 
as  follows  : 

1.  number  the  nodes  by  the  Cuthill-McKee  algorithm:  and 

2.  reorder  the  nodes  by  reversing  the  node  numbering  obtamed  m  Step  1. 

To  reverse  the  ordering  of  n  nodes  requires  only  n  operations.  Therefore,  the  overall  complexity 
of  the  RCM  algorithm  remains  bounded  by  0(m!E!>. 

It  has  been  proved  by  Liu  and  Sherman  that  the  RCM  algorithm  is  never  inferior  to  the 
Cuthill-McKee  algorithm16.  In  its  application  to  finite  element  ordering,  they  also  found  that  the 
RCM  algorithm  is  particularly  superior  when  finite  elements  of  higher  orders  are  used. 

A  listing  of  computer  subroutines  for  the  RCM  algorithm  can  be  found  in  reference  12.  A 
more  detailed  description  of  the  algorithm  can  also  be  found  in  that  reference. 
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A  TWO-STEP  APPROACH  TO  FINITE  ELEMENT  ORDERING 

In  this  section,  we  present  a  ''two-step"  approach  to  finite  element  ordering.  The  ordering 
process  is  divided  into  two  separate  tasks.  The  first  task  orders  the  finite  elements  in  the  finite 
element  mesh.  The  second  task  then  labels  the  nodal  variables. 

Representation  of  Finite  Element  Connectivity 

The  method  to  be  described  requires  an  explicit  representation  of  the  connectivity  of  the  finite 
elements  in  the  mesh.  The  connectivity  of  the  finite  elements  can  be  topologically  represented  as 
a  graph;  we  call  this  graph  a  finite  element  connectivity  graph,  or.  simply,  a  connectivity 
graph.  The  nodes  m  the  connectivity  graph  correspond  to  the  finite  elements  m  the  mesh.  The 
edges  of  the  connectivity  graph  are  used  to  describe  the  interconnections  between  the  finite 
elements.  One  possible  way  to  define  the  interconnections  between  the  finite  elements  is  to  say 
that  two  finite  elements  are  adjacent  if  they  share  a  common  node  in  the  mesh.  This  definition, 
however,  may  lead  to  a  huge  number  of  edges  in  the  connectivity  graph.  Since  the  efficiency  of 
most  existing  ordering  algorithms  is  a  function  of  the  number  of  edges  in  the  graph,  this 
representation  of  finite  element  connectivity  may  suffer  a  significant  drawback  in  terms  of  the 
execution  time  required  for  ordering  the  finite  elements.  Here,  we  examine  an  alternative 
definition  in  which  the  number  of  edges  in  the  finite  element  connectivity  graph  can  be  vastly 
reduced. 

In  finite  element  methods,  a  continuum  is  subdivided  by  imaginary  lines  or  surfaces  into  a 
number  of  finite  elements.  These  elements  are  assumed  to  be  interconnected  at  the  vertices 
situated  on  their  boundaries.  Therefore,  it  is  simple  to  assert  that  the  finite  elements  are 
topologically  interconnected  by  their  boundaries.  A  finite  element  of  n  dimensions,  where 
n  =  1.2.3,  posseses  boundaries  of  (n-l>  dimensions.  For  instance,  the  boundaries  of  a  3- 
dimensional  (volumetric)  finite  element  are  the  2~dimensional  boundary-faces,  so  that  m  a  3- 
dimensional  continuum,  the  volumetric  finite  elements  are  interconnected  by  their  boundary-faces. 
In  the  same  fashion.  2-dimensional  (planar)  finite  elements  are  bounded  and  interconnected  by 
their  i -dimensional  boundary-lines,  and  1 -dimensional  (linear)  finite  elements  are  connected  to 
their  adjacent  elements  through  the  0-dimensional  boundary-nodes.  Hence,  in  a  finite  element 


mesh,  two  finite  elements  are  said  to  be  adjacent  if  they  are  connected  at  their  boundaries,  rather 
than  at  their  vertices.  Using  this  definition,  the  nodes  in  the  finite  element  connectivity  graph  are 
the  finite  elements  in  the  finite  element  mesh,  while  the  edges  connecting  the  nodes  in  the 
connectivity  graph  correspond  to  the  imaginary  boundaries  of  the  finite  elements  separating  the 
continuum.  Examples  of  finite  element  meshes  and  their  associated  connectivity  graphs  are  shown 
in  Figure  3.  For  a  two-dimensional  (planar)  continuum,  it  is  interesting  to  note  that  the  finite 
element  connectivity  graph  is  analogous  to  the  dual  of  a  planar  graph5,  with  the  exterior  node  in 
the  dual  graph  omitted.  In  general,  the  element-element  connectivity  relationship  is  topologically 
equivalent  to  the  definition  of  a  dual  of  an  n-dimensional  complex53-  3. 

In  some  cases,  the  connectivity  of  finite  elements  cannot  be  completely  represented  using  the 
definition  given  above.  As  shown  in  Figure  4(a),  n-dimensional  finite  elements  may  not 
necessarily  be  connected  through  all  their  ( n- 1  > -dimensional  boundaries.  Another  example  may  be 
that  the  adjacent  finite  elements  do  not  have  the  same  geometric  dimensions,  as  illustrated  in 
Figure  4(b).  Nevertheless,  the  transformation  process  described  is  still  valid  and  applicable  to 
these  cases;  the  resulting  finite  element  connectivity  graph  may  at  worst  become  a  disconnected 
graph.  Each  connected  component  in  the  connectivity  graph  can  be  labelled  independently  by  the 
method  presented.  The  nodes  at  the  interface  between  two  connected  components  must  then  be 
numbered  higher  than  all  other  nodes  in  the  two  components.  The  examples  just  described  are  not 
very  common  in  practice,  and  will  not  be  pursued  further. 

There  are  several  major  advantages  in  defining  the  connectivity  of  the  finite  elements  by  means 
of  their  boundaries,  instead  of  their  nodes.  First,  the  number  of  edges  in  the  connectivity  graph 
is  minimized.  Furthermore,  this  definition  completely  disregards  the  nodal  variables,  or  vertices, 
in  the  finite  element  mesh.  Therefore,  the  ordering  of  finite  elements  can  be  performed  before 
the  types  of  finite  elements,  and  thus  the  number  of  nodal  variables,  are  chosen. 
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Ordering  the  Finite  Elements 

Once  the  adjacency  structure  of  a  connectivity  graph  has  been  established,  the  Cuthill-McK.ee 
algorithm  can  be  applied  to  number  the  nodes  of  the  connectivity  graph,  which  correspond  to  the 
finite  elements  in  the  mesh.  An  example  of  a  5  by  5  regular  mesh  with  triangular  finite 
elements  is  shown  in  Figure  5(a).  The  mesh  is  first  transformed  into  its  connectivity  graph 
shown  in  Figure  5  o).  The  nodes  in  the  connectivity  graph  are  ordered  by  the  Cuthill-McK.ee 

algorithm.  The  resulting  ordering  of  the  finite  elements  is  given  in  Figure  5(c). 

As  noted,  earlier,  the  time  complexity  of  the  Cuthill-McKee  algorithm  is  a  function  of  the 
number  of  edges  in  a  graph.  For  a  given  finite  element  mesh,  the  number  of  edges  in  the 

associated  connectivity  graph  is  considerably  less  than  the  number  of  edges  m  its  finite  element 
solution  graph  defined  previously.  When  higher  order  finite  elements  are  used,  the  difference  is 
even  more  dramatic,  since  the  number  of  edges  in  a  finite  element  solution  graph  increases 

combinaiorially  with  the  number  of  nodal  variables  per  finite  element.  On  the  other  hand,  the 
number  of  edges  in  the  finite  element  connectivity  graph  remains  constant  for  a  given 

discretization  of  the  finite  element  mesh. 

Let  the  finite  element  mesh  consist  of  ne  finite  elements,  and  let  nb  denote  the  number  of 

i 

boundaries  of  finite  element  i.  The  overall  time  complexity  of  the  Cuthill-McKee  algorithm, 

following  from  the  previous  discussion,  is  0<m!E!>.  where  m  is  the  largest  value  of  nb . 

I 

i  *  1,  ....  ne,  and  IE!  is  the  total  number  of  internal  boundaries  interconnecting  the  finite 

elements. 

Ordering  the  Nodal  Variables 

A  good  ordering  of  the  finite  elements  does  not  automatically  guarantee  a  good  node  ordering, 
in  the  sense  that  the  number  of  fill-in  entries  is  minimized.  The  ordering  of  the  finite  elements 
does  not  complete  the  ordering  of  K.  as  the  nodal  variables  associated  with  the  element  do  not 
acquire  labels  by  this  process.  Therefore,  one  must  also  consider  the  ordering  of  the  nodal 
variables.  The  proposed  strategy  is  to  label  the  nodal  variables  of  the  finite  elements  following 
the  order  in  which  the  finite  elements  are  numbered.  For  each  finite  element,  the  nodal  variables 
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are  ordered  according  to  their  valencies,  where  the  valency  of  a  node  is  the  number  of  finite 
elements  incident  at  that  node.  The  nodal  valency  is  an  indicator  of  the  degree  of  a  node. 

The  local  ordering  scheme  to  label  the  nodal  variables  is  summarized  as  follows  : 

1.  Determine  the  nodal  valency  for  each  node  in  the  finite  element  mesh. 

2.  (Mam  loop)  Enter  the  finite  elements  one  at  a  time  following  the  order  in  which  the 
finite  elements  were  previously  numbered  by  applying  the  Cuthill-McK.ee  algorithm  to 
the  connectivity  graph.  For  each  finite  element,  find  all  unlabelled  nodes  connected 
to  it  and  number  them  sequentially  in  descreasmg  order  of  their  valencies. 

3.  (Reverse  ordering)  Reverse  the  ordering  of  the  nodal  variables  obtained  in  Step  2. 

In  Step  2  of  the  local  ordering  scheme,  a  nodal  variable  with  minimum  valency  among  the 
unlabelled  nodal  variables  m  a  finite  element  is  numbered  last.  The  node  ordering  obtained  is 
then  reversed  in  Step  3.  The  motivation  behind  these  two  steps  is  that  the  nodal  variables  with 
the  lowest  valency  in  a  finite  element  will  be  eliminated  first.  This  strategy  of  minimum  degree 
(valency)  ordering  has  been  employed  in  many  popular  node  ordering  schemes,  such  as  the 
minimum  degree  ordering  algorithm1*  and  the  Cuthill-McK.ee  algorithm4,  although  the  strategy  is 
used  m  a  different  manner  m  different  schemes. 

It  has  been  mentioned  that  by  reversing  the  node  ordering  obtained  from  the  Cuthill-McKee 
algorithm,  the  result  can  be  improved  significantly.  This  property  is  also  true  m  the  local 
ordering  scheme  proposed.  First,  reversing  the  numbering  of  the  nodal  variables  which  are 
ordered  in  decreasing  order  of  their  valencies  ensures  that  internal  nodal  variables  are  numbered 
before  the  external  variables  of  the  same  element,  so  that  the  internal  nodal  variables  are 
eliminated  before  the  external  nodal  variables  of  the  same  finite  element.  Hence,  there  will  be  no 
fill-in  created  when  eliminating  the  internal  nodal  variables.  This  ordering  strategy  also  has  the 
property  that  the  nodes  situated  along  a  boundary  are  numbered  and  eliminated  before  the  corner 
nodes  of  the  same  boundary.  This  is  because  the  corner  nodes  are  m  general  connected  to  more 
finite  elements  than  the  nodes  located  along  a  boundary.  The  nodal  numbering  produced  by  the 
proposed  scheme  for  a  simple  example  of  a  5  by  5  regular  mesh  consisting  of  10-node  triangular 
finite  elements  is  shown  in  Figure  6.  The  two  properties  just  described  are  clearly  demonstrated 
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by  this  example. 

For  step  1,  given  the  element-node  incidence  table  --  a  listing  of  the  nodes  incident  on  the 

finite  elements  —  the  determination  of  the  nodal  valencies  can  be  done  in  exactly 
ne 

]  T 1  *  Z  nv  operations, 
i*l  1 

where  nv  is  the  number  of  nodal  variables  of  the  ilb  finite  element,  ne  is  the  number  of  finite 

i 

elements  in  the  mesh,  and  !  T 1  is  the  size  of  element-node  incidence  table. 

To  implement  step  2,  a  linear  insertion  may  be  used  to  sort  the  unlabelled  variables  of  each 

finite  element.  For  some  constant  c,  the  sorting  of  nv  elements  by  linear  insertion  requires 

c(nv  ->  operations1.  Thus,  the  time  complexity  of  Step  2  of  the  algorithm  requires  at  worst 
ne  ne 

c(  Z  nv:)  <  c<  Z  nv )  nv  =  c'Tjnv  operations, 

.  ,  t  —  .  ,  i  max  max 

1*1  1*1 

where  nv  *  max(nv ),  i=l,...,ne.  The  number  of  operations  required  to  reverse  the  node 

max  i 

ordering  equals  to  the  number  of  nodal  variables,  on.  in  the  finite  element  mesh.  Hence,  the 

time  complexity  for  numbering  the  nodes,  given  the  ordering  of  the  finite  elements,  is  bounded  by 

(0(  ITtnv  )  *  nn>. 

max 


EXPERIMENTAL  RESULTS 

In  this  section,  we  report  some  experimental  results  comparing  the  reverse  Cuthill-McKee 
algorithm  and  the  two-step  ordering  algorithm.  Four  finite  element  models  have  been  selected  as 

examples.  They  are  :  (1)  a  5  by  5  regular  mesh  with  triangular  finite  elements;  (2)  an  L- 

shaped  mesh  with  triangular  finite  elements;  (3)  a  cross-shaped  mesh  with  rectangular  finite 
elements:  and  (4)  a  2  by  2  by  2  model  consisting  of  rectangular  block  finite  elements.  These 
models  are  shown  in  Figure  7.  For  each  model,  linear,  quadratic  and  cubic  finite  elements  have 
been  used.  There  are  altogether  twelve  test  cases. 

The  computer  subroutines  for  the  RCM  algorithm  have  been  adopted  from  the  Sparspak  package 
developed  by  George  and  Liu  at  the  University  of  Waterloo9  and  listed  in  Reference  12.  These 

subroutines  include  finding  a  pseudo-peripheral  node  of  a  graph  and  ordering  the  nodes  by  the 


m 
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RCM  algorithm.  The  input  information  to  these  subroutines  is  the  number  of  nodal  variables  and 
an  adjacency  structure  pair  representing  the  node  connectivity  of  the  solution  graph.  A  node 
adjacency  structure  is  illustrated  in  Figure  8ia>. 

For  the  two-step  ordering  scheme,  the  same  set  of  computer  subroutines  from  Sparspak  has  been 
employed  to  number  the  finite  elements  by  the  Cuthill-McKee  algorithm,  except  that  the  step  for 
reverse  ordering  has  been  omitted.  In  this  case,  the  input  data  for  the  adjacency  structure  pair  is 
the  node  connectivity  of  the  finite  element  connectivity  graph  An  example  of  this  data  structure 
is  given  in  Figure  8< b).  Once  the  finite  elements  have  been  ordered,  the  adjacency  structure  can 
be  discarded.  The  subroutme  for  performing  the  second  step,  that  of  ordering  the  nodal  variables, 
is  listed  in  Appendix  I.  For  this  subroutine,  the  element-node  incidence  table  is  assumed  to  be 
available. 

The  test  cases  have  been  run  on  a  DECsystem  20  computer  at  Carnegie- Mellon  University.  For 
each  test  case,  the  execution  times  for  the  two-step  ordering  scheme  and  the  RCM  algorithm  are 
reported  in  Table  1.  To  minimize  timing  errors  due  to  a  multi-programmed  operating  system 
environment,  the  test  cases  have  been  run  when  the  computer  was  lightly  loaded.  As  the  results 
indicate,  the  two-step  ordering  scheme  requires  much  less  execution  lime  than  the  RCM  node 

ordering  algorithm,  except  for  the  cases  wrhen  linear  triangular  finite  elements  are  employed.  For 
those  cases,  the  execution  times  for  ordering  the  nodes  by  the  RCM  algorithm  and  for  ordering 
the  finite  elements  by  the  Cuthill-McKee  algorithm  are  approximately  the  same.  The  deficiency  of 
the  two-step  ordering  scheme  is  due  to  the  second  step  of  labelling  the  nodal  variables.  For 

cases  when  higher  order  finite  elements  are  used,  very  large  amounts  of  savings  m  execution 
times  can  be  achieved. 

For  each  test  case,  the  structures  of  the  matrix  factor  L  resulting  from  the  RCM  and  the  two- 
step  ordering  algorithms  are  summarized  m  Table  2.  in  terms  of  the  profile  of  the  stiffness 
matrix,  the  number  of  fill-m  entries  and  the  size  of  the  matrix  factor.  For  almost  all  cases,  the 
profiles  resulting  from  the  two-step  ordering  scheme  are  slightly  larger  than  those  obtained  from 
the  RCM  algorithm.  On  the  other  hand,  the  two-step  ordering  scheme,  in  most  cases,  leads  to 

less  fill-in  than  the  RCM  algorithm.  This  result  is  particularly  true  for  finite  element  meshes 
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where  higher  order  finite  elements  are  used. 

DISCUSSION 

In  this  paper,  a  "two-step"  finite  element  ordering  scheme  has  been  introduced.  In  addition  to 
the  efficiency  achieved,  the  scheme  is  also  highly  adaptable  to  various  solution  methods  commonly 
used  in  fmiie  element  analysis. 

The  major  characteristic  of  the  t«'o-step  ordering  scheme  is  that  ordering  the  finite  elements 

and  ordering  the  nodal  variables  are  separated  into  two  independent  tasks.  The  two  tasks  can  be 
implemented  as  two  separate  modules.  This  property  of  modularity  provides  flexibility  in  the 
software  design  of  finite  element  programs.  For  instance,  one  can  choose  to  number  the  finite 
elements  once  the  discretization  of  the  finite  element  mesh  is  established,  even  without  the 
knowledge  about  the  types  of  finite  elements  to  be  used.  in  view,  of  current  software 
developments  m  generating  finite  element  meshes  using  graphic  preprocesssors.  the  task  of 

numbering  the  finite  elements  can  easily  be  incorporated  as  an  additional  module  in  the  mesh 
generator  with  very  little  extra  cost. 

The  authors  recognize  that  the  two-step  ordering  scheme  presented  requires  t«o  sets  of  input 
data:  one  to  describe  the  adjacency  of  the  finite  elements  and  the  second  for  the  element-node 
incidence  table.  These  two  sets  of  data  can  both  be  generated  by  a  mesh  generator. 

Conventionally,  however,  only  the  element-node  incidence  table  if  created.  In  Appendix  11.  »e 
present  a  computer  program  to  determine  the  element-element  adjacency  structure  pair  using  the 
element-node  incidence  table  as  input.  This  computer  program  serves  primarily  for  illustration 
purposes,  and  is  not  meant  to  be  efficient.  The  authors  emphasize  that  the  finite  element 
adjacency  structure  should  be  generated  by  the  mesh  generator,  particularly  when  a  graphical 
preprocessor  is  used. 

Since  the  birth  of  the  finite  element  methods,  efforts  have  been  made  to  develop  finite  elements 
using  polynomial  interpolation  functions  of  higher  orders.  The  experimental  results  favor  strongly 
the  use  of  the  two-step  ordering  scheme  with  higher  order  finite  elements.  In  many  finite 

element  problems,  the  finite  elements  are  progressively  modified  by  using  interpolation  functions 
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of  increasingly  higher  order  so  as  to  improve  the  accuracy  of  the  solution.  At  each  modification, 
the  structure  of  the  global  stiffness  matrix  changes,  but  the  number  and  the  distribution  of  the 
finite  elements  in  the  mesh  often  remain  unchanged.  Examples  of  these  problems  are  quality 
control  m  finite  element  analysis  and  nonlinear  structural  analysis21.  For  this  type  of  problem, 
the  ordering  of  finite  elements  needs  only  be  determined  once.  The  result  can  be  used  to  reorder 
the  nodal  variables  upon  each  modification  made  to  the  finite  elements.  Moreover,  the  step  in 
ordering  the  nodal  variables  by  the  two-step  ordering  scheme  can  be  performed  much  faster  than 
a  complete  re-numbermg  by  the  node  ordering  schemes. 

In  large  scale  structural  analysis,  substructuring  is  often  used  to  divide  the  structure  into  two 
or  more  smaller  components,  called  substructures,  which  are  interconnected  at  their  boundaries. 
Topologically,  the  substructures  can  be  treated  in  the  same  manner  as  finite  elements.  Therefore, 
the  strategy  for  ordering  the  finite  elements  may  well  be  applied  to  the  ordering  of  the 
substructures.  In  substructuring  analysis,  the  external  nodal  variables  located  along  the  boundaries 
of  the  substructures  are  ordered  last.  As  discussed  previously,  this  characteristic  is  also  embedded 
in  the  node  ordering  strategy  of  the  two-step  scheme.  Hence,  there  is  no  reason  that  the  two-step 
approach  cannot  be  extended  to  an  ordering  scheme  for  substructuring  methods. 

One  of  the  characteristics  of  the  frontal  solution  method  is  that  the  assembling  of  the  finite 
elements  and  the  solution  process  are  performed  simultaneously.  For  this  solution  method,  it  is 

the  numbering  of  the  finite  elements  that  really  matters.  However,  a  proper  ordering  of  the  nodal 
variables  for  each  finite  element  can  further  reduce  the  number  of  fill-in  entries  and.  thus,  the 

computation  cost  of  the  solution  process.  In  the  two-step  scheme  presented,  the  strategy  is  to 

label  the  nodal  variables  following  the  ordering  of  the  finite  elements;  then  the  final  node 
ordering  is  reversed.  As  a  result,  the  order  in  w-hich  the  nodal  variables  are  to  be  eliminated 
ought  to  follow  the  reverse  order  m  which  the  finite  elements  are  numbered.  Therefore,  by 
entering  the  finite  elements  in  the  reverse  order  of  their  numberings,  the  assembling  and 

factorization  tasks  can  also  be  performed  simultaneously.  However,  unlike  the  frontal  solution 
method,  where  the  nodal  variables  art  arranged  in  an  arbitrary  order,  the  nodal  variables  art  pre- 
ordered  by  the  two-step  ordering  scheme  to  reduce  the  number  of  fill-in  entries.  This  pre- 


ordering  of  the  nodal  variables  has  the  advantage  that  the  data  structure  for  the  matrix  factors 
can  be  set  up  independently.  Throughout  the  entire  ordering  scheme,  the  structure  of  the  global 
stiffness  matrix  need  not  exist.  With  the  nodal  variables  pre-ordered.  one  can  proceed  directly  to 
construct  the  data  structure  for  the  matrix  factors  using  existing  symbolic  factorization  algorithms 
s  1?.  Consequently,  while  this  two-step  ordering  scheme  can  improve  the  software  modularity  for 
finite  element  programs,  many  characteristics  of  the  “natural"  frontal  solution  method  can  still  be 
retained  in  the  solution  process. 

The  authors  do  not  claim  that  the  two-step  ordering  scheme  proposed  in  this  paper  is  optimal 
in  minimizing  the  number  of  fill-in  entries.  In  fact,  it  has  recently  been  proved  that  the  problem 
of  computing  the  minimum  fill-in  is  NT-complete24.  This  ordering  scheme  is  recommended 
because  of  its  efficiency,  modularity  and  flexibility,  and  Us  potential  application  to  various  other 
finite  element  problems.  Most  of  all.  the  two-step  ordering  scheme  includes  Ihe  consideration  o! 
the  underlying  topological  structure  of  the  finite  element  mesh,  and  may  therefore  be  regarded  as 
a  "natural"  finite  element  ordering  scheme. 
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(c)  The  ordering  of  finiic  elements 


Fiyure  5:  Ordering  the  Finite  Elements  bv 
the  Cuthill- McKee  Algorithm 
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(a)  5  bv  5  regular  mesh  with  triancular  finite  elements 
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(iii)  cubic  element 
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.  Fijurc  7:  Finite  Element  Models 
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(i)  Finite  element 
solution  grapn 
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(iii)  Node-node  adjacency  structure  pair  (IA,JA) 
(a)  Adjacency  structure  pair  of  finite  element  solution  graph 
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(ii)  Element  connection  table 


Element  No.  I  II  III  XV 


M 

H 

SMI 

■■1 

m 

■ 

1 

■ 

■1 

■ 

i 

I 

JA 

II 

I 

ran 

B 

1  IV  | 

89 

(iii)  Element-element  adjacency  structure  pair  (IA,JA) 

(b)  Adjacency  structure  pair  of  finite  element  connectivity  graph 


Figure  8:  Examples  of  Adjacency  Structure  Pair 
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Appendix  I.  Fortran  Computer  Program  for  Labelling  the  Nodal  Variables 
Given  the  ordering  of  the  finite  elements,  the  subroutine  listed  below  labels  the  nodal  vanabl 
in  the  finite  element  mesh. 


C/«**«**«**************«************««***«**«*«***«**««********»*»******» 

c/«  * 

C/«  Subroutine  NODORD  :  * 

C/*  Given  the  ordering  of  the  finite  elements,  ELPM,  « 

C/*  this  subroutine  NODORD  labels  the  nodal  variables  * 

C/*  as  follows  * 

C/*  1.  Determine  the  valency  cf  each  node,  given  the 

C/*  element-r.ode  incidence  table  INC  * 

0/*  2.  label  the  nodes  cf  each  finite  element,  ir.  decreasing  * 

C/*  order  of  their  valencies  (linear  insertion  is  used)  * 

C/»  3.  Reverse  the  node  ordering  obtained  from  step  2.  * 

C/«  * 

C/«  Input  Variables  :  * 

C/«  INC  -  node-element  INCidence  table  * 

C/*  ELPM  -  Element  Permutation  vector  * 

C/«  NV  -  Number  of  Vertices  ir.  the  finite  element  mesh  * 

C/*  NVE  -  Number  cf  Vertices  per  finite  Element  « 

C/«  NE  -  Number  of  finite  Elements  in  the  mesh  * 

C/«  * 

C/*  Output  Variables  : 

C/*  PERK  -  node  PERKutatior.  vector  * 

C/*  « 

C/*  working  Variables  :  * 

C/*  IN VP  -  INVrse  Permutation  vector  * 

C/*  NDVL  -  NoDe  Valency  vector  « 

C/*  * 

C/*********************X******K**K***A****X******************************- 

SUBROUTINE  NODORD ( E1PK ,  INC,  PERK,  INV?,  NDV1,  NV,  NVE,  NE) 
IMPLICIT  INTEGER  (k-Z) 

DIMENSION  INC (NE, NVE),  PERM ( NV ) ,  INV?(NV),  EIPM(NE),  NDVl(NV) 

C/« 

C/«  Initialization 

C/» 

DO  10  K  =  1 ,  NV 
PERM (K)  =  0 
INVP(K)  =  0 
NDVL (K)  =  0 
10  CONTINUE 

0/ * 

C/*  Determine  node  valency 

C/* 

DO  20  I  =  1 , NE 
DO  20  K  =  1 , NVE 
NODE  =  INC ( I ,K) 

IF(NCDE.EQ.O)  GOTO  20 
NDVL (NODE)  =  NDVL (NODE)  +  1 
20  CONTINUE 

0/* 

C/*  Enter  the  elements  sequentially  following  their  crderir.gs 


NXTNOD  =  0 


DO  150  I  =  1,NE 
ELEM  =  ELPM(I) 

START  =  NXTNOD  +  1 
DC  110  K  =  l.NVE 
NODE  =  INC (ELEM,  K) 

IF (INVP (NODE)  .NE.O.OR.NODE.EQ.O)  GOTO  110 
NXTNOD  =  NXTNOD  -r  1 
PERM (NXTNOD)  =  NODE 
INVP (NODE)  =  NXTNOD 
110  CONTINUE 

0/* 

C/*  Local  node  ordering  :  nodes  are  sorted  in  decreasing 
C/*  order  of  their  valencies  using  linear  insertion 

C/« 

IF (START  .GE.  NXTNOD)  GOTO  15C 
K  =  START 
120  L  =  K 

K  =  K  +  1 
PERMK  =  PERM (K) 

NDVLK  =  NDVL( PERMK) 

130  IF (L  .LT.  START)  GOTO  140 

PERML  =  PERM (L) 

NDVLL  =  NDVL( PERML) 

IF (NDVLL  .GE.  NDVLK)  GOTO  140 
PERM (L+l)  =  PERML 
L  =  L  -  1 
GOTO  130 


140 

PERM (L+l) 

=  PERMK 

IF (K  .LT. 

NXTNOD) 

GOTO  120 

150 

CONTINUE 

C/« 

c/* 

Now,  reverse 

the  node 

ordering 

c/* 

KMID  =  NV/2 

J  =  NV 


DO  160  K  =  l.KMID 
I  =  PERM (K) 

PERM (K)  =  PERM ( J) 
PERM ( J)  =  I 
J  =  J  -  1 
160  CONTINUE 
RETURN 
END 
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Appendix  II.  Fortran  Computer  Program  for  Generating  the  Adjacency  Structure  Pair  of 
Element -Element  Connectivity 

Given  the  element-node  incidence  table,  the  following  set  of  subroutines  determines  the 
adjacency  structure  pair  of  element-elemen!  connectivity  of  an  arbitrary  finite  element  mesh. 
These  subroutines  are  included  only  for  illustrating  how  the  element-element  adjacency  may  be 
generated  if  the  mesh  generator  produces  only  the  element-node  incidence  table.  The  subroutines 
are  not  efficient  in  storage,  in  that  space  is  provided  for  the  entire  symmetric  node-node 
adjacency  table,  denoted  by  NEL.  Since  this  table  is  very  sparse,  the  function  HASH  couic  be 
replaced  by  one  that  exploits  sparsity1  (and  the  variable  NEDGE  reduced  accordingly >.  or  other 
efficient  techniques  for  storing  sparse  tables  may  be  used”. 


0/ 

c/> 

c/> 

c/« 

c/» 


XXXXX*XXXXXXXX*XXXXXXXXXXXXXXXXXX» 


txxxxx*xxx*x**x**xx**xxxxx»xxxxxxxxx 


Subroutine  SETUP 

Purpose  :  SETUP  is  the  main  subroutine  tc  generate  the 
element-element  adjacency  structure  from  the 


0/* 

element-node  incidence  table. 

INO,  for  a  planar 

X 

c/« 

finite  element  mesh. 

X 

c/* 

X 

c/« 

Input  variables  : 

X 

c/» 

INC  -  element-node  incidence  table 

X 

0/* 

X 

c/« 

Global  Variables  : 

* 

c/* 

NE  -  no.  of  elements  in  the  mesh 

X 

c/« 

NV  -  no.  cf  nodes  in  the  mesh 

X 

c/« 

NVE  -  no.  of  nodes  oer  element 

X 

o/« 

NEPLS1  -  (NE  -  1) 

* 

c/* 

NEDGE  -  possible  maximum  r.o.  cf  edge 

s  ir.  the  solution  graph 

X 

0/ * 

’(  =  NV  «  (NV-1)  /  2  ) 

X 

c/* 

EDGKAM  -  total  no.  cf  edges  ir.  the  s 

clutior.  graph 

* 

0/ * 

EDGFL1  -  (EDGMAX  +  1) 

X 

c/» 

NELS  -  counter  for  the  edge-elemer.t 

structure  pair 

* 

0/* 

MAXSU3  -  max.  subscript  used  in  the 

element -element 

* 

c/» 

adjacency  structure  pair  (see  subroutine  ELMAOJ) 

* 

c/» 

X 

c/« 

Output  : 

X 

0/ * 

The  element-element  adjacent  structu 

re  pair  (IA,CA)  is  stereo 

* 

0/. 

in  the  array  PCCL  from  location  LOCI A 

to  ( LOO J  A- MAMS  U  5 - 1 )  . 

* 

0/. 

The  array  IA  is  stored  in  PCOL(LOCIA) 

to  POOL (LOCIA-NE) ; 

X 

c/* 

the  array  JA  is  stored  ir.  PCCL(LOCUA) 

to  PCCL l ICC JA-MAXSUE-1) . 

X 

C/*  * 

C/»  Subroutines  used  :  * 

C/»  EDGSET,  EDG5LT ,  ELMADJ  * 

C/« 

£/***»****************************************«************************** 

SUBROUTINE  SETUP ( INC ,  POOL) 

IMPLICIT  INTEGER  (A-Z) 

COMMON  /SYSTEM/  N'E ,  NV ,  NVE,  NEPLS1 
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C/» 

C/« 

c/< 

c/» 


c/« 
c/» 
c/« 
c/ * 


COMMON  /GRAPH/  NEDGE,  EDGMAX,  EOGPLi ,  NELS ,  MAXSUH 
COMMON  / ICSET/  IN,  10,  DEBUS 
EGG I CAL  DEBUG 

DIMENSION  INC ( NE , NVE) ,  ?00L(1) 

Phase  1  :  to  set  up  the  edge-element  adjacency  structure 
LOCNEL  =  1 

LOCSTA  =  LOCNEL  +  NEDGE 

CALL  EDGSET ( INC,  POOL (LOCNEL) ,  POOL (LOCSTA) ) 

Phase  2  :  to  build  the  edge-element  adjacency 
structure  pair 

EDGPLl  =  EDGMAX  *  1 
LOCEL  =  LOCSTA  +  EDGPLl 
LCCPOS  =  LOCEL  -  NELS 

CALL  EDGBLT(INC,  POOL (LOCNEL) ,  POOL (LOCSTA) ,  PCOL(LCCEL), 
POOL (LOCPOS ) ) 

Phase  3  :  to  build  the  element-element,  adjacency 
structure  pair 


LOCI A  =  LCCPOS 

LOCO A  =  LOCI A  +  NEPLS1 

CALL  ELMADJ (INC,  POOL (LOCNEL) ,  POOL ( LOCSTA) ,  PCCL(LOCEL), 

*  POOL ( LOCIA) ,  POOL (LGCJA) ) 

RETURN 

END 

C/ »»***«»*****«******««***''****«*«*******»***•***«*«***»»*«»«>"<*******«* 

C/« 

C/«  Function  HASH 

C/*  Purpose  :  to  determine  the  location  of  an  edge  in 

C/«  the  symmetric  node  adjacency  matrix 

C/* 

C/« 

c/« 
c/» 
c/« 
c/» 
c/* 
c/» 

INTEGER  FUNCTION  HASH(N1,N2) 

IMPLICIT  INTEGER  (A-Z) 

J  =  AMAXC ( N1 , N2 ) 

I  =  AMINO ( NI , N2 ) 

HASH  =  (  ( J-2 ) * ( J-l) ) /2  -  I 

RETURN 

END 

C/ *«»*«**************»*»*«»»***««******«*««****«********««»*»«***«*»*»«*• 
C/*  * 

C/*  Subroutine  EDGSET  « 

Purpose  :  to  set  up  counter  START  of  the  edge-element  * 

adjacency  structure  pair  (START,  EL)  * 

{EL  will  be  determined  in  subroutine  EDGELT . }  * 


Input  Variables  : 

Nl,  N2  -  nodes  in  the  mesh  (or  graph) 

Output  Variable  : 

HASH  -  location  of  edge  (N1,N2)  in  the  symmetric  node 
adjacency  matrix 


c/* 

c/« 

c/* 

c/« 

c/« 


Inout  Variable  : 
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C/«  INC  -  element-node  incidence  cacie 

C/« 

C/»  Output  Variables  : 

C/ *  NEL (HASH (Nl, h'2) )  -  edge  no.  assigned  to  node  pair  (N1,N2) 

C/«  {  Only  edges  with  more  than  one  incident  element 

C/»  are  assigned  a  number.  } 

C/*  START  -  starting  position  for  the  edge-element  adjacency 

C/*  structure  pair 

C/»  {  START()c+l)  -  STARTOc)  -  no.  of  elements  adjacent 

C/»  to  edge  X  } 

C/» 

C/*  Intermediate  Variable  : 

C/»  NEL ( HASH ( Nl , N  2 ) )  -  in  step  1,  the  array  is  a  node-ad ;acency 

C/»  matrix  denoting  number  of  elements  into  edge  (N1.N2; 

c/* 

c/ ***«*****«***»********«****«**«***»««»*****«*****«**««**«•»»******« 

SU3RCUTINE  EDGSET ( INC,  NEL,  START) 

IMPLICIT  INTEGER  (A-Z) 

COMMON  /SYSTEM/  NE,  NV,  NVE,  NEPLS1 

COMMON  /GRAPH/  HEDGE,  SDGMAX,  EDGPL1 ,  NELS ,  MAXSl'3 

DIMENSION  I NO (NE, NVE) ,  NEL (NEDGE) ,  START (1) 

C/* 

C/«  Initialization 

C/» 

DC  5  I  =  1, NEDGE 
NEL (I)  =  0 
5  CONTINUE 

C/* 

C/«  Step  1  :  to  build  the  node  adjacency  matrix 

C/« 

DO  1C  I  =  1, NE 
NN  =  NVE 

DO  10  J  =  1,  NN-1 
Nl  =  INC (I , J) 

IF(Nl.EQ.C)  GOTO  10 
DC  10  K  =  J+l , NN 
N2  =  INC ( I , K) 

IF  (Nl .  EQ . N2  .OR.  N2.EQ.0)  GOTO  10 
EDGE  =  HASH (Nl , N2 ) 

NEL (EDGE)  =  NEL (EDGE)  *  1 
10  CONTINUE 

C/» 

C/*  Step  2  : 

C/*  (1)  to  set  up  counter  START  for  the 

C/*  edge-element  adjacency  structure 

C/*  (2)  to  assign  a  number  to  the  edges  (with  two 

C/*  or  more  incident  elements  ) . 

c/» 

POSIT  =  0 
NELS  =  I 

DC  15  I  =  1,  NEDGE 

IF (NEL (I)  .LE.  1)  GOTO  14 
POSIT  =  POSIT  +  1 
START (POSIT)  =  NELS 
NELS  =  NELS  +  NEL (I) 

NEL ( I )  =  POSIT 


NEL ( I )  = 


14 


0 


IS  CONTINUE 

EDGMAX  =  POSIT 

START (EDGMAX+i)  =  NELS 

RETURN- 

END 

:/**»»**««• . . . . 

c/« 

C/»  Subroutine  EDGELT 

C/«  Purpose  :  to  build  the  array  ED  of  the  edge-element  adjacency 

C/*  structure  pair  (START,  EL). 

c/« 

0/*  Input  Variables  : 

C/«  INC  -  element-node  incidence  table 

0/*  NEL(k)  -  edge  no.  assigned  to  node  pair  (Ni,N2) 

0/*  {  k  =  HASH  ( N1 ,  N2 )  } 

c/*  START  -  starting  position  for  the  edge-element  adjacency 

C/*  structure  pair 

0/*  {  START ( k+1 )  -  START ( k)  =  no.  of  elements  adjacent 

C/«  to  edge  k  } 

c/« 

C/«  Output  variable  : 

C/*  ED  -  edge-element  adjacency  structure  pair 

C/*  {  The  set  of  elements  incident  to  edge  it  is  stored 

C/»  in  ED  (START  (It)  )  to  ED  (START  (k*l) )  .  ) 

C/» 

C/»  working  Variable  : 

0/*  POSIT (k)  -  pointer  to  current  empty  position  in  EL  for  edge  k 

c/» 

c/»  ********************************************************************  * 
SUBROUTINE  EDGBLT (INC,  NEL,  START,  EL,  POSIT) 

IMPLICIT  INTEGER  (A-Z) 

COMMON  /SYSTEM/  NE,  NV,  NVE ,  NEPLS1 
COMMON  /GRAPH/  NEDGE,  EDGMAX,  EDGPL1 ,  NELS,  MAXSUB 
DIMENSION  INC (NE , NVE) ,  NED (NEDGE) ,  START (EDGPL1 ) ,  EL (NELS) , 

POSIT (EDGMAX) 

Initialize  working  vector  POSIT 

DO  10  I  =  1,  EDGMAX 
POSIT (I)  =  START (I) 

CONTINUE 

Fill  edge-element  adjacency  EL 

DO  2C  I  =  1,NE 
NN  =  NVE 

DC  20  J  *  1,  NN-1 
N1  =  INC ( I , J) 

IF(Nl.EQ.O)  GOTO  20 
DC  20  K  =  J+l,  NN 
N2  =  INC ( I ,K) 

IF (N2 . EQ . 0  .OR.  N1.EQ.N2)  GOTO  20 
EDGE  =  HASH (N1 ,N2 ) 

IF (NEL (EDGE)  .EQ.  0)  GOTO  20 
EDGACT  =  NED (EDGE) 

EL  ( POS I T  ( EDGACT )  )  =  I 
POSIT (EDGACT)  =  POSIT (EDGACT)  *  1 


Inode  Nl! 


inode  N2I 
: check  self  locp: 
I  location  m  NED  I 

ledge  number! 
lassigr.  I  to  ED: 

I  update  PCS I Tier: 


C/* 

c/* 

c/» 


c/» 

0/* 

c/* 


20  CONTINUE 
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RETURN 

END 

C/********************************************************************** 

C/» 

C/«  Subroutine  ELMADJ 

C/«  Purpose  :  to  set  up  the  element-element  adjacency  structure 

C/«  pair  ( IADJ, JADJ) . 

c/« 

C/*  Input  Variables  : 

C/«  INC  -  element-node  incidence  table 

C/»  NELOO  -  edge  no.  assigned  to  node  pair  (N1,N2) 

C/*  {  k  =  HASH (N1,N2)  } 

C/*  START  -  starting  position  for  the  edge-element  adjacency 

C/«  structure  pair 

C/*  {  START (k+1)  -  START (k)  =  no.  of  elements  adjacent 

C/*  to  edge  k  ) 

C/*  EL  -  edge-element  adjacency  structure  pair 

C/«  {  The  set  of  elements  incident  to  edge  k  is  stored 

C/«  in  EL (START (k) )  to  EL (START (k+1) ) .  } 

C/« 

C/»  Output  variables  : 

C/*  (IADJ, JADJ)  -  element-element  adjacency  structure  pair 

C/« 

c/********************************************************************** 

SUBROUTINE  ELMADJ (INC,  NHL,  START,  EL,  IADJ,  JADJ) 

IMPLICIT  INTEGER  (A-Z) 

COMMON  /SYSTEM/  NE,  NV,  NVE,  NEPLS1 
COMMON  /GRAPH/  NEDGE,  EDGMAX,  EDGPL1 ,  NSLS,  MAXSU3 
COMMON  /IOSET/  IN,  IC,  DEBUG 
LOGICAL  DEBUG 

DIMENSION  INC  (NE,  NVE)  ,  NEL  (NEDGE)  ,  START  (EDC-PL1)  ,  EL(NELS), 

*  IADJ  (NEPLS1)  ,  .fADJ(l) 

C/* 

C/»  Initialize  IADJ(1) 

C/* 

IADJ (1)  =  1 

C/» 

C/*  For  each  element,  find  all  incident  elements  and 

C/*  aueue  them  in  JADJ 


DO  50  I  =  1 , NE 
SLMENT  =  I 
LI  =  IADJ (I) 

FILL  =  0 
NN  =  NVE 

DO  40  J  =  1,  NN-1 
NI  =  INC ( I , J) 

IF(Ni.EQ.O)  GOTO  40 
DC  4C  K  =  J+l,  NN 
N2  =  INC ( I , K) 

IF (Nl . EQ. N2  .OR.  N2.EQ.0)  GOTO  40 
EDGE  =  HASH (Nl ,N2) 

IF (NEL (EDGE) .EQ.C)  GOTO  40 
EDGN’JK  =  NEL  (EDGE) 

BEGIN  =  START (EDGNUK) 

END  =  START (EDGNUM-1)  -  1 
DC  30  IEL  =  BEGIN,  END 
ELM  =  EL (IEL) 


!for  element  I! 


Inode  Nl! 


!node  N2! 

! check  self  loop! 
’.location  in  NEL! 

ledge  number! 

!for  all  elements 
!  incident  to 
:  edge  edgnum: 

!a  neigh,  element 
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IF ( ELM . £Q . ELMENT )  GOTO  30 

L2  =  LI  +  FILL 

IF(FILL.EQ.O)  GOTO  20 

!same  as  element  I 

DO  10  JJ  =  LI,  L2-1 

check  if  there 

IF (JADJ (JJ) . EQ.ELM)  GOTO  30 

is  a  duplication 

10 

CONTINUE 

in  lest  : 

20 

JADJ (L2)  =  ELM 

FILL  =  FILL  +  1 

queue  it  m  JADJ! 

30 

CONTINUE 

40 

CONTINUE 

IADJ (ELMENT+1)  =  I AD J (ELMENT)  +  FILL 

.'update  IADJ! 

50 

C/» 

CONTINUE 

MAXSUB  =  IADJ (NEPL31) 

IF( .N0T.DE3UG)  RETURN 

C/» 

c/« 

Print  out  results 

WRITE (10, 1100) 

DO  50  I  =  1 ,  NE 

LI  =  IADJ  (I) 

L2  =  IADJ  ( I-fl)  -  1 

IF(L2.LT.L1)  GOTO  60 

WRITE (10, 10C0)  I,  LI,  L2 ,  (JADJ(K),  K=L1,L2) 

50 

CONTINUE 

RETURN 

1000 

FORMAT (T5 ,  14,  TIC,  16,  T18,  16,  (T30,i0Io)) 

1100 

FORMAT (1H1,  T2 ,  'ELEMENT',  T10,  ’ IDAJ ( I ) ' ,  TIB, 

IADJ ( I+l) ' , 

* 

T30 ,  'JADJ  (NEIGHBORING  ELEMENTS) 

END 
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